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Abstract 

We consider the implementation of the split-step method where the linear 
part of the nonlinear Schrodinger equation is solved using a finite-difference dis- 
cretization of the spatial derivative. The von Neumann analysis predicts that 
this method is unconditionally stable on the background of a constant-amplitude 
plane wave. However, simulations show that the method can become unstable 
on the background of a soliton. We present an analysis explaining this instabil- 
ity. Both this analysis and the instability itself are substantially different from 
those of the Fourier split-step method, which computes the spatial derivative by 
spectral discretization. We also found that the modes responsible for the numer- 
ical instability are supported by the sides of the soliton, in contrast to unstable 
modes of linearized nonlinear wave equations, which (the modes) are supported 
by the soliton's core. 
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1 Introduction 



The split-step method (SSM), also known as the operator-splitting method, is widely 
used in numerical simulations of evolutionary equations that arise in diverse areas of 
science: nonlinear waves, including nonlinear optics and Bose-Einstein condensation 
PQ [5], atomic physics [6J [T] , studies of advection-reaction-diffusion equations [8] [TT]. 
and relativistic quantum mechanics [12]. In this paper we focus on the SSM applied 
to the nonlinear Schrodinger equation (NLS) 

in* — (3u xx + ju\u\ 2 = 0, u(x, 0) = uq(x). (1) 

Although the real-valued constants (3 and 7 in (CQ) can be scaled out of the equation, 
we will keep them in order to distinguish the contributions of the dispersive (u xx ) and 
nonlinear (m|m| 2 ) terms. Without loss of generality we will consider 7 > 0; then solitons 
of CQ) exist for < 0. 

The idea of the SSM is that (Tj[|) can be easily solved analytically when either the 
dispersive or the nonlinear term is set to zero. This alows one to seek an approxi- 
mate numerical solution of (CQ) as a sequence of steps which alternatively account for 
dispersion and nonlinearity: 



for n from 1 to n max do: 

u(x) = u n (x) exp (i-f\u n (x)\ 2 At) (nonlinear step) 

, solution of iu t = /3u xx at t = At . 
u n+ i(x) = < (dispersive step) 

with initial condition u(x, 0) = u(x) 



(2) 



where the implementation of the dispersive step will be discussed below. In ([2]), At 
is the time step of the numerical integration and u n (x) = u(x,nAt). Scheme ([2]) 
can yield a numerical solution of ([I]) whose accuracy is O(At). Higher-order schemes, 
yielding more accurate solutions (e.g., with accuracy 0(At 2 ), 0(At 4 ), etc.), are known 
[T3| [HI E], but here we will restrict our attention to the lowest-order scheme ([2]); see 
also the paragraph after Eq. fl33|) below. 

The implementation of the dispersive step in (T2]) depends on the numerical method 
by which the spatial derivative is computed. In most applications, it is computed by 
the Fourier spectral method: 

u n+1 (x) = I' 1 [exp(i/3k 2 At) F[u{x)} ] . (3) 

Here J 7 and J 7 ™ 1 are the discrete Fourier transform and its inverse, k is the discrete 
wavenumber: 

- 71/ Ax < k < tc/Ax, (4) 

and Ax is the mesh size in x. However, the spatial derivative in (T5]) can also be com- 
puted by a finite-difference (as opposed to spectral) method [151 EE E3 EH] • For exam- 
ple, using the central-difference discretization of u xx and the Crank-Nicolson method, 
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the dispersion step yields: 

. < + i~u m / C + + 1 1 -2< +1 + <y- 1 1 - 2u m + u 

% At 2 V Ax 2 Ax 2 

where = u(x m , nAt), x m is a point in the discretized spatial domain: — L/2 < x m < 
L/2, and L is the the length of the domain. Recently, solving the dispersive step of (j2J) 
by a finite-difference method has found an application in the electronic post-processing 
of the optical signal in fiber telecommunications [19]. That post-processing must be 
done in real time, and hence the speed of its implementation becomes a key factor. 
Ref. [20] showed that finite-difference implementations (referred to there as the Finite 
and Infinite Impulse Response filters; see its Eqs. (13) and (22)) of the dispersive step 
of (jSJ provide considerable speed improvement over the Fourier-based solution (j3J); see 
also Sec. 2.4.2 in [3]. Finally, let us note that the version of the NLS where the second 
derivative is replaced by its finite-difference approximation, as on the right-hand side 
(r.h.s.) of ([S]), is of considerable interest also in its own right |21j . 
In what follows we assume periodic boundary conditions: 

u{-L/2,t) = u{L/2,t); (6) 

the case of other types of boundary conditions is considered in Appendix A. We will 
refer to the SSM with spectral flH]) and finite-difference (jHJ) implementations of the 
dispersive step in (J2]) as s-SSM and fd-SSM, respectively. Our focus in this paper will 
on the fd-SSM. 

Weideman and Herbst [15] used the von Neumann analysis to show that both 
versions, s- and fd-, of the SSM can become unstable when the background solution of 
the NLS is a plane wave: 

u pw = (A/yfy) e luJ ^\ A = const, oo pw = \A\ 2 . (7) 

Specifically, they linearized the SSM equations on the background of (j7|): 

Un = Upw + U nj \u n \ < \u n \ (8) 

and sought the numerical error in the form 

u n = Ae xtn ~ ikx , A = const. (9) 

The SSM is said to be unstable when for a certain wavenumber k one has: (i) Re(A) > 
in (|9"1), but (ii) the corresponding Fourier mode in the original equation ([1]) is linearly 
stable. Weideman and Herbst found that the s- and fd-SSMs on the background (|7j) 
become unstable when the step size At exceeds: 

At thr , s w Ax 2 /(ti\(3\), for s-SSM © & © (10) 
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and 

At thrifd = Ax/^2\/3\ \A\ 2 only for (3 > , for fd-SSM (J2) & © (11) 

respectively. Note that for f3 < 0, the fd-SSM simulating a solution close to the plane 
wave (E]) is unconditionally stable. Typical dependences of the instability growth rate, 
Re(A) > 0, on the wavenumber is shown in Fig. [TJ Let us emphasize that the SSM is 
unstable for At > At thr even though both its constituent steps, ([2]) and either (J3]) or 
(EJ), are numerically stable for any At. 




Figure 1: Growth rate of numerical instability of the s-SSM (a) and fd-SSM (b) on the 
plane-wave background. The dotted horizontal line indicates how the maximum growth 
rate depends on the wave's amplitude. In (a), k m7T , m = 1, 2, . . . are the wavenumbers 
where the mth resonance condition holds (see [21]): l/^/c^At = ititc. 

Solutions of the NLS (and of other evolution equations) that are of practical interest 
are considerably more complicated than a const ant- amplitude solution (J7J). To analyze 
stability of a numerical method that is being used to simulate a spatially varying 
solution, one often employs the so-called "principle of frozen coefficients" [22J (see also, 
e.g., [231 EG])- According to that principle, one assumes some constant value for the 
solution u and then linearizes the equations of the numerical method to determine the 
evolution of the numerical error (see (IE]) and (|9])). However, as we show below, this 
principle applied to the SSM fails to predict, even qualitatively, important features of 
the numerical instability. 

As a first step towards developing a general framework of stability analysis on the 
background of a spatially varying solution, we analyzed [23] the instability of the s-SSM 
on the background of a soliton of the NLS: 

u sol {x, t) = A^2h sech (Ax/y/=P) e ioj «° lt = U sol (x) e l ^\ u sol = A 2 . (12) 

First, we demonstrated numerically that the instability growth rate in this case is very 
sensitive to the time step At and the length L of the spatial domain; also, its depen- 
dence on the wavenumber is quite different from that shown in Fig.[]Ja). Moreover, the 
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instability on the background of, say, two well-separated (and hence non-interacting) 
solitons can be completely different from that on the background of one of these solitons. 
To our knowledge, such features of the numerical instability had not been reported for 
other numerical methods. Moreover, they could not be predicted based on the principle 
of frozen coefficients. We then demonstrated that all those features could be explained 
by analyzing an equation satisfied by the numerical error of the s-SSM: 



where v(x, t) is proportional to the continuous counterpart of u n (x) = u(x, nAt) defined 
in (jSJ), and k n is defined in the caption to Fig. [IJ Note that ffTB"]) is similar, but not 
equivalent, to the NLS linearized about the soliton: 



The extra fc^-term in ([TBI indicates that the potentially unstable numerical error of 
the s-SSM has a wavenumber close to k n . 

In this paper we present an analysis of the numerical instability of the fd- (as op- 
posed to s-) SSM on the soliton background. This analysis is based on an equation for 
the numerical error which, as ({TBI , is a modified form of the linearized NLS. However, 
both that equation and its analysis are qualitatively different from those [23] for the 
s-SSM. Moreover, the modes that render the s- and fd-SSMs unstable are also qual- 
itatively different. Namely, for the s-SSM, the numerically unstable modes contain 
just a few Fourier harmonics and hence are not spatially localized. On the contrary, 
the modes making the fd-SSM unstable are localized and are supported by the sides 
(i.e., tails) of the soliton. To our knowledge, such "tail-supported" localized modes, as 
opposed to those supported by the soliton's core, have not been reported before. 

The main part of this manuscript is organized as follows. In Sec. II we present 
simulation results showing the development of numerical instability of the fd-SSM 
applied to a soliton. In Sec. Ill we derive an equation (a counterpart of (fTBl) ) governing 
the evolution of the numerical error, and in Sec. IV obtain its localized solutions that 
grow exponentially in time. In Sec. V we summarize the results. In Appendix A we 
show how our analysis can be modified for boundary conditions other than periodic. In 
Appendix B we describe the numerical method used to obtain the localized solutions 
in Sec. IV. In Appendix C we discuss how the instability sets in. 

2 Numerics of fd-SSM with soliton background 

We numerically simulated Eq. (TjQ) with j3 = — 1, 7 = 2, and the periodic boundary 
conditions (jHJ) via the fd-SSM algorith (j2J) & (JHD- The initial condition was the soliton 
(TT2]) with A=l: 



iv t - u sol v - f3(v xx + k\v) + 7|w so i| 2 (2v + v*) = 0, 



(13) 



iut - u soi u - (3u xx + 7|m so1 | 2 (2m + u*) = 0. 



(14) 



Uq(x) = sech (x) +£(#); 



(15) 
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the noise component £(x) with zero mean and the standard deviation 10~ 10 was added 
in order to reveal the unstable Fourier components sooner than if they had developed 
from a round-off error. 

Below we report results for two values of the spatial mesh size Ax = L/N, where 
N is the number of grid points: Ax = 40/2 9 and Ax = 40/2 10 . We verified that, for 
a fixed Ax, the results are insensitive to the domain's length L (unlike they are for 
the s-SSM [24] ) as long as L is sufficiently large. Also, at least within the range of 
Ax values considered, the results depend on Ax monotonically (again, unlike for the 
s-SSM). 

First, let us remind the reader that the analysis of [15] on a constant-amplitude 
background ([ZD for f5 < predicted that the fd-SSM should be stable for any 
For the soliton initial condition ( ITBl and the parameters stated above, our simulations 
showed that the numerical solution becomes unstable for At > Ax. For future use we 
introduce a notation: 

C = (At/ Ax) 2 . (16) 

In Fig.^a) we show the Fourier spectrum of the numerical solution of (CEJ), (TT51) obtained 
by the fd-SSM with C = 1.05 (i.e., slightly above the instability threshold) at t = 1400. 
The numerically unstable modes are seen around the end points of the spectral axis. 
At t — 1400, these modes are still small enough so as not to cause visible damage 
to the soliton: see the solid curve in Fig. [2](b) . However, at a later time, the soliton 
begins to drift: see the dashed line in Fig. Mjo), that shows the numerical solution at 
t = 1800. Such a drift may persist over a long time: e.g., for C = 1.05, the soliton still 
keeps on moving at t ~ 4000. However, eventually it gets overcome by noise and loses 
its identity. 




Figure 2: See explanation in the text. 

1 Recall that the stability of instability of the numerical method is in no way related to that of the 
actual solution. In fact, the plane wave (0 is well-known to be modulationally unstable for j3 < 0, 
while it is modulationally stable for f3 > 0. 
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We observed the same scenario for several different values of Ax, L, and C (for 
C > 1). The direction of the soliton's drift appears to be determined by the initial 
noise; this direction is not affected by the placement of the initial soliton inside the 
spatial domain. The time when the drift's onset becomes visible decreases, and the 
drift's velocity increases, as C increases. 

The soliton's drift is a nonlinear stage of the development of the numerical instabil- 
ity. It will be explained in Sec. V. In the linear stage, the numerically unstable modes 
are still small enough so that they do not visibly affect the soliton or one another. To 
describe this stage, we computed a numeric approximation to the instability growth 
rate Re(A) defined in (jUJ): 

^ / max | | for \ / noise floor 

I k ~ k max at time= t J \ at time= 

Re (A) | computed = " , (17) 

where fc max = n/Ax (see PJ). The so computed values of the instability growth rate 
are shown in Fig. [3] along with the results of a semi-analytical calculation presented in 
Sec. IV. 




C-1 

Figure 3: Growth rate of the numerical instability for Ax = 40/2 9 (solid line — analysis 
of Sec. IV, stars — computed by (fT7l) ) and for Ax = 40/2 10 (dashed line — analysis 
of Sec. IV, circles — computed by ifPTj) ). 

The above numerical results motivate the following three questions: (i) explain the 
observed instability threshold At (see the sentence before ( fl6l) ); (ii) identify the modes 
responsible for the numerical instability; and (iii) calculate the instability growth rate 
(see Fig. [3]). In Sec. IV we will give an approximate analytical answer to question (i). 
However, answers to questions (ii) and (iii) will be obtained only semi-analytically, i.e., 
via numerical solution of a certain eigenvalue problem. 
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3 Derivation of equation for numerical error for fd- 
SSM 



Here we will derive a modified linearized NLS that is a counterpart of (1131) . which 
was obtained for the s-SSM in [24]. In view of periodic boundary conditions (0), the 
finite-difference implementation (EJ of the dispersive step in (j2J) can be written as 



u n+1 (x) 



T- 1 [e iP{k) 7 [u{x)\ ] , (U 



,- P / M 1 + 2z(3rsm 2 (kAx2) , ,„ . . N1 At 

e lP(fc) = ^ = exp 2iarctan(2/3r sin 2 (k Ax 2) ) , r = — — , 

1 -2i/3r sin 2 (fc Ax/2) F[ v v 7 7 ;J ' Ax 2 ' 

(19) 

where J 7 , J 7-1 were defined after Q. For |fcAx| 1, the exponent in (Tl9|) equals that 
in (J3]). However, for |/cAx| > 1, they differ substantially: see Fig. 0J It is this difference 
that leads to the instabilities of the s- and fd-SSMs being qualitatively different. 



/ s-SSM 




fd-SSM 



Figure 4: Normalized phase: \(3\k 2 At for the s-SSM (dashed) and as given by ( 1191) for 
the fd-SSM (solid). In both cases, r = 5. The horizontal line indicates the condition 
of the first resonance: |-P(/c)| = 7r. 



Using Eqs. ([2]) and (Tl8|) . one can write, similarly to Eq. (3.1) in [24], a linear 
equation satisfied by a small numerical error u n of the fd-SSM: 



+i 



3 iP(*) 



J 7 



3 i 7 KI 2 At/,~ 



(«n + «7At(«b M n + kb| 2 Mn) ) 



(20) 



Here M n is defined similarly to (jSJ), with w b being either w pw or u so \, depending on the 
background solution. The exponential growth of u n can occur only if there is sufficiently 
strong coupling between u n and u* in f |20|) . This coupling is the strongest when the 
temporal rate of change of the realtive phase between these two terms is the smallest. 
In [24] we showed that this rate can be small only for those k where the exponent P(k) 
is close to a multiple of it. Using ( TT9l) (see also Fig. H]), we see that this can occur only 
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for sufficiently high k where sin 2 (/cAx/2) = 0(f) rather than 0(Ax 2 ). Then: 

P ^ k) = {3rsm 2 (kAx/2) + ° 

f (k - k max ) 2 Ax 2 ^( 1 ((k-k max )AxY , . 
" * fir 4/3r + U \r* + r 1 ' 

where fc max = tt/Ax. We have also used that 

r = At/ Ax 2 = C/At > f , (22) 

given that the numerical instability was observed in Sec. 11 for C = 0(1). 

In order for the truncation of the Taylor expansions used in (|21|) to be self-consistent, 
one needs to assume that 

0(l/r 2 ) < 0{(k - k max ) 2 Ax 2 ) < 0(1); 

note that 0(l/r 2 ) = 0(Ax 2 ). We arbitrarily declare (k — fc max ) 2 Ax 2 to be in the 
middle of that range: (k — /c max ) 2 Ax 2 = O(Ax), which implies that the range of 
wavenumbers for which the instability develops satisfies: 

|fc-fcmax| = 0(l/\^) < k max = n/Ax. (23) 

This turns out to be consistent with the result shown in Fig. 12(a). 

Substituting the first three terms on the r.h.s. of ( l2"Tj) into (]2"U|). using f l2"2~j) . and 
introducing a new variable 



(e™) n u n = (-l) n u n) (24) 



one obtains: 



T[v n+1 ) = exp -— <^ f + ^ -!■ \ x 
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nKI 2 A* |~ n + ilAt{ ulv* n + \u h \ 2 v n )} . (25) 



Note that (T251) describes a small change of v n occurring over the step At, because for 
At — > 0, the r.h.s. of that equation reduces to J 7 ^]- Therefore we can approximate 
the difference equation ([215]) by a differential equation. To that end, first recall from 
( |23l) that the wavenumbers of v n are on the order of fc max ; hence we seek 



v n (x) = e lfemaxX w n (x). (26) 

The effective wavenumber of w n is then (k — fc max ), and according to ff23|) w n varies 
slowly over the scale O(Ax). Introducing the scaled space and wavenumber by 

X = x/e, k = (k - fc max )e, e = Ax/2, (27) 
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one rewrites (|25|) as: 

T sc [w n+ i] = exp (-^{1 + « 2 }J ^"sc [e i7|Mb|2A * {w n + iiAt(ulw* + \u h \ 2 w n )} 

(28) 

where now J-" sc is the Fourier transform with respect to the scaled variables ( 1271) . In 
handling the v* term in (125]) . we have used the fact that on the spatial grid x m = mAx, 
one has: 

v*(x m ) = e~ ik ^ Xm w* n {x m ) = e- i7Tm w* n (x m ) = e inm w* n (x m ) = e ik ™* Xm w* n (x m ). 



Next, the s-SSM (J2]), ([3]) can be written as 



(29) 



When |/3 1 k 2 At < 1 and ~/\u\ 2 At <C 1, this is equivalent to the NLS ([1]) p/ws a term 
proportional to 

At [pdl, 7|M| 2 ]„M + 0(At 2 ), (30) 

where [. . . , . . .]_ denotes a commutator (see, e.g., Sec. 2.4 in [3]). Equation (|28|) has 
the form of a linearized Eq. (129]) with a different coefficient in the dispersion term and 
with an extra phase. Therefore, (1281) must be equivalent to a modified linearized NLS, 
with the modification affecting only the corresponding terms: 

iw t + (w xx - w)/(C/3) + 7(«b** + 2K| 2 w) = 0, (31) 



plus a term proportional to the linearized form of the commutator ( 1301) . Neglecting 
that latter term as small (of order O(At)) compared to the rest of the expression and 
denoting ip — w exp(— iu^t), we rewrite (|3~T|) as: 

iiPt + 6i> + *P xx /{C(1) + 7 t/ b 2 (e X ) (2V + r) = 0, (32) 

where 

<f=-Wb-l/(C73). (33) 

Here Wb is either a; pw or w so i, and is either constant or U so \, depending on whether 
the background solution is a plane wave ([7]) or a soliton (TT2|) . The modified linearized 
NLS (1321 for the fd-SSM is the counterpart of Eq. ([TBI that was derived for the s-SSM. 

Our subsequent analysis of the instability of the first-order accurate fd-SSM (T2]) & 
(J5j) will be based on Eq. ( 1321) . The instability of the second-order accurate version of 
this method, where the order of the nonlinear and dispersive steps is alternated in any 
two consecutive full time steps [13J, is the same as that of the first-order version. The 
instability of higher-order versions (e.g., 0(At 4 )-accurate) can be studied similarly to 
how that was done in Ref. [241 for the s-SSM. 
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The boundary conditions satisfied by if) are still periodic: 

if)(-L/(2e),t) = if;(L/(2e) 7 t). (34) 

This follows from the fact that u n (x) satisfies the periodic boundary conditions and 
from ( 126|) . given that for k max = it/ Ax and L/2 = MAx with some integer M, 

There are three differences between Eq. (J32l) and the linearized NLS fTT4l) . Most 
importantly, (l32il has the opposite sign of the dispersion term. This is explained by 
the shape of the curve P(k) for the fd-SSM in Fig. H]at high wavenumbers, where the 
curvature is opposite to that at ~ 0. Secondly, unlike the (— w so i)-term in (1141 . the 
5-term in fl32|) with < can be either positive or negative, depending on the value 
of C. Thirdly, the "potential" £7^ (ex) (when U\, = U so i) is a slow function of the scaled 
variable x- That is, solutions of ( |32|) that vary on the scale x — 0(1) "see" the soliton 
as being very wide. This should also be contrasted with the situation for the s-SSM, 
where the modes described by Eq. (fl~3|) "see" the soliton as being very narrow [24J. 

Before proceeding to find unstable modes of Eq. (1321) with U\> = U so \, let us note 
that (1321) with Ub = const confirms the result of Ref. [33] regarding the instability of the 
fd-SSM on the plane-wave background. Namely, for (3 < 0, Eq. fl32l) with = const 
describes the evolution of a small perturbation to the plane wave in the modulationally 
stable case (see, e.g., Sec. 5.1 in [3]). That is, for (3 < 0, there is no numerical instability, 
in agreement with [15] . On the other hand, for /3 > 0, Eq. ( 1321) describes the evolution 
of a small perturbation in the modulationally unstable case, and hence the plane wave 
of the NLS ([1]) can become numerically unstable. The corresponding instability growth 
rate found from (132!) and Eq. (5.1.8) of [3] can be shown to agree with the one that 
can be obtained from Eq. (37) and the next two unnumbered relations in [15] . An 
example of this growth rate is shown in Fig. [T](b). Also, using our (1321) and Eq. (5.1.8) 
of Ref. [3], the threshold value of At can be shown to be given by (TIT]) , in agreement 
with [15]. 

4 Unstable modes of modified linearized NLS for 
fd-SSM 

In this section we focus on the case where the background solution is a soliton; hence 
f3 < and = U so y (see f[T2|) ). Substituting into fl32l) and its complex conjugate the 
standard ansatz |25j (?/>(x,£), V ; *(X)^)) — (0l(x)> 02(x)) eA< an d using yet another 
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rescaling: 



A _ 2A x „ Cf3\_ 2 { 1 



A = ^A, V(y) = 2Cp 2 sech 2 (y), 



(35) 



d 2 x + D - V{eX) \\\\\4= * W, ( 36 ) 



one obtains: 



where 03 = diag(l, —1) is a Pauli matrix, = (0i, 2 ) , and T stands for a transpose. 
If (0, A) is an eigenpair of (EBJ), then so are (<Ji<j>, -A), (0*, -A*), and (<Ti<j>*, A*), 
where 

/ 1 
1 




0i 



is another Pauli matrix. Note also that A is defined in the same way as in (jH]); hence 
Re(A) ^ indicates an instability. Below we will use shorthand notations A_r = Re(A) 
and A 7 = Im(A). 

We begin analysis of ( )36l) with two remarks. First, this equation is qualitatively 
different from an analogous equation that arises in studies of stability of both bright 
[25] and dark [26] NLS solitons in that the relative sign of the first and third terms 
of f )36|) is opposite of that in [251 126]. This fact is the main reason why the unstable 
modes supported by ( l36l) are qualitatively different from unstable modes of linearized 
NLS-type equations, as we will see below. While the latter modes are supported by 
the soliton's core (see, e.g., Fig. 3 in [27]), the unstable modes of ( 1361) are supported 
by the soliton's "tails" . 

Second, from ( 136]) and (135]) one can easily establish the minimum value of parameter 
C where an instability (i.e., A# ^ 0) can occur. The matrix operators on both sides 
of ( )36l) are Hermitian; the operator o~3 on the r.h.s. is not sign definite. Then the 
eigenvalues A are guaranteed to be purely imaginary when the operator on the left- 
hand side (l.h.s.) is sign definite [28]; otherwise they may be complex. The third term 
on the l.h.s. of ( 136]) is negative definite, and so is the first term in view of (|34p . The 
second term, D, is negative when 

C < 1/(|/3|A 2 ). (37) 

Thus, ( ITB|) and fl37j) yield the stability condition of the fd-SSM on the background of 
a soliton. We will show later that an unstable mode indeed first arises when C just 
slightly exceeds the r.h.s. of ( |37l) . 

Since the potential term in (13"6"]) is a slow function of X, it may seem natural to 
employ the Wentzel-Kramers-Brillouin (WKB) method to analyze it. Below we show 
that, unfortunately, the WKB method fails to yield an analytic form of unstable modes 
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of (156]) . However, it still indicates where they can exist. Away from "turning points" 
(see below) the WKB-type solution of (156]) is: 



= ( a+ e 0+/t + b + e- e+/e ) 0+ + (a_e e - /e + 6_e^- /e ) <p_, (3* 
where a±, 6± are some constants, and 



>2 



-D + 2^±v / l /2 -A 2 , ^ = F(eX), e' = d9/d(eX) 7 (39a) 



VA ± VA 2 - V< 

<P± = " 777 I , I . (39b) 

' - py 2 (K 2 -A 2 )] 1/4 ^ -VAtVA 2 -\/ 2 

At a turning point, say, X = X , the solution ( 138|) . (159]) breaks down, which can 
occur because the denominator in ( 139bj) vanishes. In such a case, one needs to obtain 



a solution of (156]) in a transition region around the turning point by expanding the 
potential: V(eX) = V(eX ) + e(X — X ) V'(eX ) + . . ., and then solving the resulting 
approximate equation. For a single Schrodinger equation, a well-known solution of 
this type is given by the Airy function. This solution is used to "connect" the so far 
arbitrary constants a±, b± in ( 158]) on both sides of the turning point. 

Now a turning point of ( 156]) is where: either (i) 6' + = or 9'_ = 0, or (ii) (V(eX)) 2 — 
A 2 = 0. The former case can be shown (see, e.g., [29]) to reduce to a single Schrodinger 
equation case, where the solution in the transition region is given by the Airy function. 
However, at present, no such transitional solution is analytically available in case (iij 2 ) 
[30l [31] . Therefore, the solutions (|38j) canot be "connected" by an analytic formula 
across such a turning point, and hence one cannot find the eigenvalues A in this case. 

Based on the past experience with unstable linear modes of nonlinear waves, it is 
reasonable to assume that that unstable solutions of (156]) must be localized. This allows 
one to predict the location of these solutions, which will then be verified numerically. 
Indeed, expressions (J58J) and ( 159aD show that (6' ± ) 2 > for (2V - D) > and V > A 
(i.e., inside the soliton's core). The solution ( 158]) cannot exponentially grow from, 
say, the left side of the soliton to the right. Indeed, on the scale of Eq. ( 156]) . the 
soliton is very wide and hence, if the solution (!58j) is of order one on the left, it would 
have become exponentially large on the right side of the soliton. Therefore, inside the 
soliton, the solution (J58J must be exponentially decaying. From (I39aj) it follows that 
it is also possible for a solution with A R ^ to decay outside the soliton (i.e., where 
V ~ 0) when D > 0. Thus, one can expect that a localized mode of (156]) is to be 
lumped somewhere at the soliton's side and decay in both directions away from that 
location. 

In Fig. |5]we show the first (i.e., corresponding to the greatest A^j) such a mode for 
L = 40, N = 2 9 points (hence e = Ax/2 w 0.04), A = 1, p = -1, 7 = 2. For these 



2 Note that in this case, (V 2 — A 2 ) 1 / 4 <p + and (V 2 — A 2 ) 1 / 4 ^- are linearly dependent. 
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parameters, the threshold given by the r.h.s. of (|37|) is C = 1, and parameter D in 
(136]) is related to C by: 

D = C-1. (40) 

The numerical method of finding these modes is described in Appendix B, and the 
modes found by this method are shown in Fig. GD^a) for different values of C. In 
Fig. Efb) we show the same modes obtained from the numerical solution of the NLS 
(CQ) by the fd-SSM. These modes were extracted from the numerical solution by a high- 
pass filter, and then the highest-frequency harmonic was factored out as per (1261) . The 
agreement between Figs. EJ^a) and E(b) is seen to be good. In Fig. [5](d) we show the 
location of the peak of the first unstable mode, computed both from fl36l) and from 
the numerical solution of ([1]), versus parameter C. The corresponding values of the 
instability growth rate A were shown earlier in Fig. [31 Let us stress that A for the 
localized modes of (1361) was found to be purely real up to the computer's round-off 
error (~ 10~ 15 ). There also exist unstable modes with complex A, but such modes are 
not localized and have smaller growth rates than the localized modes. 




x x x C- 1 



Figure 5: (a): Profiles of the first localized mode on the right side of the soliton for 
different values of C, as found by the numerical method of Appendix B. (b): Same 
as in (a), but found from the numerical solution of (CQ), as explained in the text, (c): 
The modes at both sides of the soliton found from the numerical solution of ([I]). Note 
that these modes do not "see" each other because of the barrier created by the soliton, 
and hence in general have different amplitudes as they develop from independent noise 
seeds. In panels (a)-(c), the potential is sech 2 (eA) (see ( |35l) ) and the amplitude of 
the mode is normalized to that of the potential, (d): Location of the peak of the first 
localized mode, found by the method of Appendix B (solid line) and from the solution 
of ([1]) (stars). Similar data for L = 40 and iV = 2 10 are very close to those in (d) and 
hence are not shown. 

As C increases from the critical value given by (1371) . the localized unstable mode 
becomes narrower and also moves toward the center of the soliton. Moreover, higher- 
order localized modes of ( l36|) arise. Typical profiles of the second and third modes 
are shown in Fig. El along with the parameter C for which such modes first become 
localized within the spatial domain. In Appendix C we demonstrate that the process 
of "birth" of an eigenmode that eventually (i.e., with the increase of C) becomes 
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localized, is rather complicated. In particular, it is difficult to pinpoint the exact value 
of parameter C where such a mode appears. Therefore, the C values shown in Fig. [6] 
are accurate only up to the second decimal place. 




x x c - 1 

Figure 6: Similar to Fig. |5^a), but for the second (a) and third (b) localized modes, (c): 
C values where localized modes of increasing order appear. Stars — for e = 40/1024, 
circles — for e = 40/2048. 



5 Conclusions 

The main contribution of this work is the extension the (in) stability analysis of the 
fd-SSM beyond the von Neumann analysis in order to include spatially-varying back- 
ground solutions, such as the soliton (fT2"j) . We showed that, as previously for the 
s-SSM [21] . this is done via a modified equation, ( )32|) . derived for the Fourier modes 
that approximately satisfy the resonance condition 

\P\k 2 At = 7i. (41) 

Note that such modified equations — (fi~3l) for the s-SSM and ( 1321) for the fd-SSM — are 
different. Their analyses are also qualitatively different, and so are the modes that are 
found to cause the instability of these two numerical methods. For the s-SSM, these 
modes are almost monochromatic (i.e., non-localized) waves ~ exp(±ikx) that "pass" 
through the soliton very quickly. It is this scattering of those waves on the soliton that 
was shown [21] to lead to their instability. In contrast, for the fd-SSM considered in this 
work, the dominant unstable modes are stationary relative to the soliton. Moreover, 
they are localized at the sides, as opposed to the core, of the soliton. To our knowledge, 
such localized modes were not reported before in studies of (in) stability of nonlinear 
waves. We consider this finding another significant contribution of this work. 

Let us stress that, as previously for the s-SSM, the principle of frozen coefficients 
fails to even qualitatively describe the instability of the fd-SSM on the soliton back- 
ground. Indeed, that principle assumes that the spatially varying coefficient U£(ex) in 
( |32l) can be approximated by a constant. However, such an approach predicts [15] that 
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any solution of the NLS with f3 < 0, including the soliton, should be stable (see (TTTj) 
and the end of Sec. Ill), contrary to the findings of this work. In fact, the unstable 
modes of (1361) localized at the soliton's sides owe their existence to the spatial variation 
of the soliton's profile. 

It was straightforward to obtain an approximate threshold, ( 1371) (where C is given 
by f|T6|) ). beyond which the numerical instability may occur. Our simulations showed 
that the instability does indeed occur just slightly above that threshold. However, 
finding the exact threshold, C cr , remains an open problem, and so is the calculation of 
the growth rate of the localized unstable modes; see Appendix C. 

Let us note that the instability of the fd-SSM was studied in a recent paper [IB"] . 
However, the focus of that paper was on a rigorous mathematical proof — by a com- 
pletely different method than we used here — of the stability of the numerical scheme 
under a certain condition, rather than on the details of the development of the in- 
stability, as in our work. That condition — their Eq. (2.9) — is a constraint on the 
parameter r = At/ Ax 2 , whereas our condition f )37|) is a constraint on the parameter 
At/ Ax. Clearly, as Ax — > 0, our constraint allows for a greater At than that of [18] . 
and hence is sharper. Also, it is consistent with the numerical experiments reported in 
that paper. 

Finally, let us show how our results can qualitatively explain the observed dynamics 
of the numerically unstable soliton — see the text after Eq. ( fT6l) . Let u unst be the 
field of the unstable modes on the soliton's sides. At an early stage of the instabilty, 
|w U nst| *^ A, the amplitude of the soliton. Also, its characteristic wavenumbers are 
much greater than those of the soliton: see Fig. [2(a) and f )23|) . Accordingly, substitute 
u = u soi + M unst into the NLS ([I]) and discard all the high-wavenumber terms to obtain: 

i( U so\)t - 0( U 8o\)xx + 7«sol|Wsol| 2 = -27«solKnst| 2 - (42) 

This is the equation for a perturbed soliton with the perturbation being, in general, 
not symmetric about the soliton's center (see Fig. [5(c)), because the modes on the 
left and right sides of the soliton do not "see" each other and hence can have different 
amplitudes. Such a perturbation is known (see, e.g., [3], Sec. 5.4.1) to cause the soliton 
to move, which is precisely the effect we observed. 
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Appendix A: Modified linearized NLS for fd-SSM 
with non-periodic boundary conditions 

We consider the Dirichlet boundary conditions (b.c.) and without loss of generality 
(for the purposes of this derivation) set them equal to zero. Neumann or mixed b.c. 
can be treated similarly, and lead to similar results. 

The equation for the dispersive step of the SSM is still given by ([SJ). However, now 
instead of (J6]) we assume: v? n+l = 0, u n+ \M = 0, where m = and m = M are the 
end points of the spatial grid. Then (J5J) can be rewritten as [32J 

(X + (0r/2)A) u n+1 = (1 - (i0r/2)A) u, (43) 

where: u = [u 1 , u 2 , . . . , u 1 ] T , similarly for u n +i, I is an (M — 1) x (M — 1) identity 
matrix, and A is an (M — 1) x (M — 1) tridiagonal matrix with (—2) on the main 
diagonal and (+1) on the sub- and super-diagonals. 

The starting point of our derivation in Sec. Ill, Eq. ffl8l) . has exacly the same form 
for the case of the Dirichlet b.c, except that J 7 is replaced with T - an expansion 
over the complete set of the eigenvectors of A; similarly, J 7-1 is replaced by T _1 . The 
exponential in ( fl9|) that acts on the jth eigenvector is replaced by 

1 ~ »/?r A 3 -/2 

I + 1^/2' 1 ] 

where Aj is the corresponding eigenvalue [32J: 

Xj = -4 sin 2 (nj/(2M)). (45) 



Equations (1441) . (1451) and the middle expression in (1191) coincide provided that we 
identify: 

k = jir/(MAx) = jrt/L. (46) 

However, we are still a step away from proving that the modified linearized NLS for 
the Dirichlet b.c. case is the same as that equation for periodic b.c. This is because 
— k 2 , which is the Fourier symbol of the second derivative, is not the symbol of the 
second derivative under the transformations T and T _1 . Under those transformation, 
the required symbol is given by (|45p . We will now use this observation to supply the 
last step and show that the modified linearized NLS for the case of Dirichlet b.c. is 
indeed the same as (13T1) . This follows from (j4"4"]) - (l4"6"l) and a calculation that is similar 
to (HQ: 

e iP « w - | 1 + ' 



i/3r(l-sin 2 ((A;-A; max )Ax/2)) 



x , 1 , sin 2 ((£;-fc max )Ax/2)) , ^ 

i/3r ifir 
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where we have used that sin(fc max Aa;/2) = 1 and that for highly oscillatory eigenvectors 
of A, one has (k — /c max )Ax <C 1. The last term on the r.h.s. of ( H7|) is the desired 
symbol of the second derivative, and then the rest of the derivation is the same as that 
leading to (1281) . From it one obtains the same modified linearized NLS as ( 13TT) . Our 
numerical simulations of the NLS using the fd-SSM with zero Dirichlet b.c. confirm 
this conclusion. 



Appendix B: Numerical solution of eigenproblem ( 1361 ) 

We work with (l36|) written in an equivalent form: 



(7 3 [d 2 x + D - V{eX) 



2 1 
1 2 



iA <J 3 = i(A - A )a 3 (f), 



(48) 



where the reason to include a constant Ao will be explained later. We discretize (1481) 
using Numerov's method. It approximates the equation &xx = F(Q,X) by a finite- 
difference scheme 

AX 2 



m+l 



Tfl I 1 



2$ m + $ 



12 



(49) 



with accuracy 0(AX 4 ); here $ m = $(X m ), F m = F($ m ,X m ), etc., and m = 
1, . . . , M — 1. Then, for the discretized solution f fc = [(fil, . . . , 0fc f_1 ] T = 1,2) 
one obtains: 



(-1) 



k-l 



AX 



;A peI + Af per {DI - 2V - iA {-l) h - 1 X} 



ffc — A/" per Vf3_fc 



= «(A-A )A/' p erffc, 

(50) 

where X is defined after f T4"5j) : A pcr is as in f T4"5j) except that its (1,M — l)th and 
(M — 1, l)th entries equal 1 (to account for the periodic boundary conditions); Af pcr 
has a similar structure as A per : (A/p e r)m,m = 10/12 (see (|49j) ). (A/per)(m-i),m = 
(A/" per ) m ,( m -i) = 1/12, (A/" P cr)i,(M-i) = (A^er) = V 12 , and the rest of its en- 
tries are zero; and V = diag(V l , . . . , V^ 1 ^ 1 ). Next, defining the combined vector and 
matrices: 



fi 
f 2 



A 



per 



"4-per O 

O A pC x 



i A/p Cr 



O Afr 



per 



V 



2V V 
V 2V 



o- 3 



O -1 



where O is the (M — 1) x (A<f — 1) zero matrix, one rewrites (!50|) as: 
1 



0"3 



-Aper + D A/per ~ A/ p0 rV - 2A CT 3 .A/per 



f = z(A - A )A/p er f . (51) 
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This equation has the form of the generalized eigenvalue problem Qi = A"Hf where % = 
A/'per is a positive definite matrix. This problem can be solved by Matlab's command 
eigs. As its options, we specified that we were looking for 24 smallest-magnitude 
eigenvalues and the corresponding eigenmodes. Beyond the instability threshold there 
are several modes with complex A. We visually inspected them and found that the 
most unstable mode was also the most localized and also had a real eigenvalue. 

We verified that the eigenvalues did not change to five significant figures whether 
we used AX = 1/10 or 1/20; so we used AX = 1/10. Finally, the need to shift the 
eigenvalues by a Ao occurred for relatively large D: D > 0.1. There, the localized 
modes were no longer among those with the 24 smallest-magnitude eigenvalues (see 
above), and significantly increasing the number of modes sought caused eigs to fail. 
Therefore, we had to shift the eigenvalues by Ao in order to make the mode of interest 
appear among the 24 smallest-eigenvalue ones. 



Appendix C: "Birth" of localized unstable mode 

The instability growth rates plotted in Fig. [3] are monotonic functions of the the pa- 
rameter C. This, however, occurs only when C is sufficiently beyond a critical value, 
C CI , where the dominant (i.e., with the greatest- |A^|) unstable mode is created. Near 
C cr , which is slightly above the threshold value given by the r.g.s. of ( |37|) . the evolution 
of the greatest- 1 A/j | eigenvalue is quite irregular. 

Below we present results about this evolution for the dominant unstable mode 
(shown in Fig. [5]) for L = 40, N = 2 9 (i.e., e = 40/1024 « 0.04) and the rest of the 
parameters being the same as listed in Sec. II, i.e.: /3 = —1, 7 = 2, and A — 1. Then, 
the control parameter C is related to the free term in (1361) by ( 1401) . While we have 
been unable to rigorously establish an analytical expression for C cr , we will present 
a hypothesis as to what it may be. The main message that we intend to convey is 
that the "birth" of a localized eigenmode of Eq. ( 1361) occurs via a complex sequence 
of bifurcations, in contrast to a single bifurcation that typically takes place when an 
unstable mode of a nonlinear wave is "born" (see, e.g., |33j). 

We numerically observed that eigenmodes of (1361) with A# ^ emerge not only 
from the origin (i.e. when imaginary eigenvalues "collide" at A = and give rise to 
two real ones), but also from "collision" of eigenvalues with Aj ^ whereby twc§ 
complex eigenvalues are "born". However, the very first (i.e., for the smallest C) 
unstable mode does emerge out of the origin. We will now show that this mode is 
essentially non-localized and, moreover, it is not the mode that eventually becomes 
the dominant unstable mode, whose growth rate is plotted in Fig. [3j The reason why 

3 As per the remark after (|36|) . there is also a pair of eigenvalues with (—A/), so in this case there 
exists a quadruplet of complex eigenvalues. 
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we still chose to discuss the former mode while being primarily interested in the latter 
one, will become clear as we proceed. 

For a mode with A = 0, Eq. ( 1361) can be split into two uncoupled Schrddinger 
equations: 

( d\ + D - V± V(eX))<j>± = 0, </>± = </>i± 02, (52a) 
where v_ = 1 and u + = 3. Note that <p± satisfy the periodic boundary conditions: 

± (-L/(2e)) = 0±(L/(2e)). (52b) 

In Fig. [7Ja) we show an example of a nontrivial solution of (|52|) . In view of the periodic 
boundary conditions, this figure is equivalent to Fig. [TJb). Recall from Sec. IV that 
the eigenmode is exponentially small inside the soliton. Then the solution shown in 
Fig. UJb) can be thought of as being localized inside the valley bounded by the two 
"halves" of the potential. Using this observation, one can estimate the isolated values 
of D for which one of the equations ( 152al) . along with (I52bl) . has a nontrivial solution, 
by the WKB method. The condition for the existence of a mode localized inside the 
valley of Fig. [7(b) is given by the Bohr-Sommerfeld formula: 

( / 1Cft + [ /( ) V / D-uV(eX)dX = tt (n+^Y (53) 

\J-L/(2e) JX Tight J V <V 

where v is either V- or u + , n is an integer, and Xi e ft )r i g ht are the turning points (see 
Sec. IV), where 

D-vV(eX Mt , light )=0. (54) 

The number of full oscillation periods of the mode inside the valley equals n; for 
example, in Fig. [7J n = 3. 




Figure 7: (Color online) (a): A solution of f )52|) (solid); potential sech 2 (eX) (red 
dotted). The amplitude of the solution is normalized to that of the potential, (b): 
Same as (a), but that panel is "cut" along the vertical dotted line at the center, and 
the resulting halves are interchanged. 
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When D < 1, the sech potential in (|54|) can be approximated by an exponential: 
sech 2 (eX) fa 4exp(— 2eX). Then, using (|35|) and fl54"l) . we reduce (1531) to 



bright 

with Xieft = — X ris ht and 



/•i/(2e) , / ,\ 

v^D / ^1 - exp[-2e(X - X right )] dX = - in + -\, 



(55) 



bright = — In — . (56) 

Neglecting the exponentially small terms of the order 0( exp[— (L — 2eX rig ht)] ), one 
obtains from ( j55|) : 

^ L _ ln ^^__2(i-l n 2)^ = e7r L+iY (57) 

Note that the WKB condition ( 1531) . and hence ( |57|) . is valid when n is sufficiently 
large. In particular, it is not supposed to accurately predict the "birth" of the first 
unstable mode, where n = 0. Indeed, Eq. (|57|) predicts that such a mode (for z/ = 1) 
emerges at D fa 5.9 ■ 10~ 6 , while numerically (see Appendix B) it is found at D fa 
1.6 • 10~ 5 . (A similar mode for v = 3 emerges at a slightly higher value of D.) Formula 
fl57|) becomes accurate to the fourth significant figure in D for n > 20. 

As we noted above, the first unstable mode is not the one that eventually becomes 
the dominant unstable mode. It disappears already at D fa 1.7 • 10~ 5 , and there is an 
adjacent interval of D values where all the eigenvalues of ( 1361) are purely imaginary 
(i.e., the soliton is numerically stable). As D increases, higher-order "real" (i.e., with 
Aj = 0) modes appear and disappear in a similar fashion, as do quadruplets of modes 
with complex A. In both these types of modes, Ar is fairly small: |A#| < D/10. There 
also exist intervals of D, of increasingly small length, where all A's are purely imaginary. 
This situation persists until the dominant unstable mode appears at D CT (= C CT — 1). 
This occurs as follows. 

First, at D fa 0.012134, a "real" mode appears (see Fig. [Ha,b)), and from this 
point on there always exists a "real" mode, even though the particular mode "born" 
at D fa 0.012134 disappears later on. Specifically, at D fa 0.012928, another "real" 
mode acquires A# greater than that of the mode "born" at D fa 0.012134, and the 
latter mode soon disappears (Fig. [8](c,d)). A similar switchover between "real" modes 
occurs at least one more time near D fa 0.013750 (not shown). Next, another "real" 
mode is "born" via a cascade of bifurcations near D fa 0.0162 (Fig. El(e,f)), and its 
A/? crosses that of the previously dominant-A^ "real" mode near D fa 0.01635. At 
D = 0.0170, these two dominant "real" modes have A R fa 1.4 ■ 10~ 3 and 1.5 • 10~ 3 
(Fig. E](e)). Finally, these two modes gradually approach each other while crossing at 
least once more near D = 0.01725. At D = 0.023 and beyond, their eigenvalues are 
the same to five significant figures. Thus, remarkably, the dominant unstable mode 
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x 10~ 4 x 10" 4 x 10~ 4 




Figure 8: (Color online) Real and imaginary parts of selected modes of fl36|) . including 
the most unstable mode, (a) & (b): 0.01200 < D < 0.01220; (c) & (d): 0.01285 < 
D < 0.01300; (e) & (f): 0.01610 < D < 0.01700. Same line colors, styles, and widths 
are used to indicate the same modes within one pair of panels (e.g., (c) & (d)). The 
same line colors/styles/ widths in different pairs of panels (e.g., in (a) & (b) and (c) & 
(d)) do not imply the same modes. 



eventually becomes doubly degenerate. We verified that the same also holds for the 
higher-order localized unstable modes. 

To conclude, we present a hypothesis as to why the value C CT , where a "real" mode 
appears permanently (see above), is near C = 1.012 — 1.013. Let us interpret ( IBTl) in 
a way that the n on its r.h.s. is not necessarily an integer, but a continuous function 
of the parameter D. For those values of D when n is an integer, a mode with a real 
A either appears or disappears at the origin A = 0. Evaluating n at the values of D 
listed in the previous paragraph in connection with Figs. E^a)-(d), one finds: 

at D = 0.012134 : n\ u=l « 29.02, n\ v= i - u\ v=3 « 0.99; (58a) 
at D = 0.012928 : n\ v=x w 30.02, n\ v=1 - u\ u=3 « 1.02; (58b) 
at D = 0.013750 : n\ v=1 w 31.04, n\ v=1 - v\ v=3 « 1.05. (58c) 

That is, both n\ v= i and n\ u=3 are simultaneously very close to integers. At D = 
0.012928, one of the "real" modes has not yet disappeared while the next one has 
appeared (Fig. [8]^c)). From (|58|) we observe that at this value of D, the difference 
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(ri|„=i — n\ u=3 ) exceeds 1 for the first time. Thus, we hypothesize that C cr = 1 + D cr is 
found from the condition that (n\ u=1 — n\ u=3 ) exceeds 1 for the first time. Verification 
of this hypothesis requires a deeper analytical insight than we have at the moment. 
Moreover, finding a value of C past which the dominant real eigenvalue increases 
monotonically (as seen in Fig. E|e)) is also an open question. 
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List of captions 

Fig. 1: Growth rate of numerical instability of the s-SSM (a) and fd-SSM (b) on the 
plane-wave background. The dotted horizontal line indicates how the maximum growth 
rate depends on the wave's amplitude. In (a), k m7T , m = 1, 2, . . . are the wavenumbers 
where the mth resonance condition holds (see [21]): l/^A^At = mix. 

Fig. 2: See explanation in the text. 

Fig. 3: Growth rate of the numerical instability for Ax = 40/2 9 (solid line - 
analysis of Sec. IV, stars — computed by fflTj) ) and for Ax = 40/2 10 (dashed line - 
analysis of Sec. IV, circles — computed by ffTTl)). 

Fig. 4: Normalized phase: \(3\k 2 At for the s-SSM (dashed) and as given by (fUl 
for the fd-SSM (solid). In both cases, r = 5. The horizontal line indicates the condition 
of the first resonance: |-P(&)| = ir. 

Fig. 5: (a): Profiles of the first localized mode on the right side of the soliton for 
different values of C, as found by the numerical method of Appendix B. (b): Same 
as in (a), but found from the numerical solution of ([I]), as explained in the text, (c): 
The modes at both sides of the soliton found from the numerical solution of ([1]). Note 
that these modes do not "see" each other because of the barrier created by the soliton, 
and hence in general have different amplitudes as they develop from independent noise 
seeds. In panels (a)-(c), the potential is sech 2 (eX) (see fl3"5l) ) and the amplitude of 
the mode is normalized to that of the potential, (d): Location of the peak of the first 
localized mode, found by the method of Appendix B (solid line) and from the solution 
of (CQ) (stars). Similar data for L = 40 and N = 2 10 are very close to those in (d) and 
hence are not shown. 

Fig. 6: Similar to Fig. |5^a), but for the second (a) and third (b) localized modes. 

(c) : C values where localized modes of increasing order appear. Stars — for e = 
40/1024, circles — for e = 40/2048. 

Fig. 7: (Color online) (a): A solution of ( !52l (solid); potential sech 2 (eX) (red 
dotted). The amplitude of the solution is normalized to that of the potential, (b): 
Same as (a), but that panel is "cut" along the vertical dotted line at the center, and 
the resulting halves are interchanged. 

Fig. 8: (Color online) Real and imaginary parts of selected modes of (136]) . including 
the most unstable mode, (a) & (b): 0.01200 < D < 0.01220; (c) & (d): 0.01285 < 
D < 0.01300; (e) & (f): 0.01610 < D < 0.01700. Same line colors, styles, and widths 
are used to indicate the same modes within one pair of panels (e.g., (c) & (d)). The 
same line colors/styles/ widths in different pairs of panels (e.g., in (a) & (b) and (c) & 

(d) ) do not imply the same modes. 
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